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We present a semiclassical calculation of the generalized form factor K a b{r) which characterizes 
the fluctuations of matrix elements of the operators a and b in the eigenbasis of the Hamiltonian of a 
chaotic system. Our approach is based on some recently developed techniques for the spectral form 
factor of systems with hyperbolic and ergodic underlying classical dynamics and / = 2 degrees of 
' freedom, that allow us to go beyond the diagonal approximation. First we extend these techniques 

to systems with / > 2. Then we use these results to calculate K a b(T~). We show that the dependence 
on the rescaled time r (time in units of the Heisenberg time) is universal for both the spectral and 
the generalized form factor. Furthermore, we derive a relation between K a b(T) and the classical 
<**] time-correlation function of the Weyl symbols of a and b. 
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Q 1 I. INTRODUCTION 
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A. Overview 



In a number of recent works the quantum spectral statistics of closed chaotic systems was investigated in the 
'— semiclassical limit. According to a conjecture by Bohigas, Giannoni and Schmit p] (BGS), the fluctuations of the 
energy levels are system-independent and coincide with the predictions of random-matrix theory (RMT) if the system 
has a chaotic underlying classical dynamics. Numerical and experimental investigations carried out on a great variety 
| of systems support the BGS conjecture 0, ■ 
r- — i ■ In the semiclassical limit, Gutzwiller's trace formula Q provides a suitable starting point for the calculation of 
spectral correlation functions. It relates the quantum mechanical density of states to a sum over classical periodic 
orbits which are characterized by an amplitude and a phase that is obtained from the action of the orbit. A prominent 
example of a quantum correlation function is the two-point energy-energy correlation function and its Fourier trans- 
form, the spectral form factor K(t). In this case a semiclassical analysis faces the problem of evaluating a double sum 
over periodic orbits which requires an appropriate quantitative treatment of classical action correlations. Averaging 
the form factor of a given system over an energy window that is large compared to the mean level spacing implies 
that only pairs of orbits with a small action difference of the order of h yield significant contributions. The leading 
contribution given by the terms with vanishing action difference is obtained within the diagonal approximation 
. This approximation accounts for all orbit pairs where an orbit is associated to itself or to its time-reversed partner if 
time-reversal symmetry is present. 

Only recently a method for a systematic inclusion of certain orbit pairs with small but nonzero action differences was 
developed for systems with time-reversal symmetry 00] ■ This approach, originally formulated in the configuration 
space, was first applied to the next-to-leading order correction of the spectral form factor of a uniformly hyperbolic 
system showing agreement with the RMT predictions. In subsequent works, extensions to a phase-space formulation 
applicable to nonuniformly hyperbolic systems H, |j| and to higher-order corrections 0] were proposed. However, 
all these previous considerations were restricted to systems with / = 2 degrees of freedom, e.g., two-dimensional 
billiards, while the more general RMT conjecture is independent of /. In this work, we present a theory which applies 
to hyperbolic and ergodic Hamiltonian systems with an arbitrary number / of degrees of freedom (/ > 2). These 
systems are characterized by a set of system-specific time scales, namely the (/ — 1) positive Lyapunov exponents. 
However, we will prove that going beyond the diagonal approximation the final result for the spectral form factor is 
independent of / and coincides in a universal way with the RMT predictions. 
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From an experimental point of view it is desirable to furthermore develop a theory that describes not only statistical 
properties of the energy spectrum but also of quantum mechanical matrix elements, as entering, for example, in cross 
sections. The fluctuations of the diagonal matrix elements of the operators a and b in the eigenbasis of the Hamiltonian 
can be described by a generalized form factor K ab {T). 

Similarly to the spectral form factor, one expects that K a b{r) also shows universal features as % — > and depends 
only on averaged classical quantities, like the averages over the constant-energy surface and the time correlation 
function of the Weyl symbols of a and b. Analytical results exclusively based on the diagonal approximation have to 
some extent confirmed this statement [lj, US, LLJ- Gil ■ ^ n this work, we generalize these results beyond the diagonal 
approximation. 

In the following two subsections we briefly recall the semiclassical theory for the form factor based on Gutzwiller's 
trace formula. In section [HJ we discuss the diagonal approximation and the origin of the off-diagonal corrections 
in the special case of the spectral form factor. We then show how the results for two-dimensional systems can be 
extended to higher-dimensional ones, leading once more to universality and agreement with RMT predictions in the 
semiclassical limit. The matrix-element fluctuations described by the generalized form factor K a b(T) are then studied 
in Section IlIII The leading-order (diagonal approximation) and next-to-leading-order terms are determined. 



B. Generalized form factor: definitions and main results 



We introduce the weighted density of states 

d a (E) =tr( a 8(E-H)) = ^(n|a|n) <5(£ - E„) (1) 



for a quantum observable a. This is in generalization of the spectral density of states where a is given by the identity 
operator, i.e., a = 1. In Eq. IjTJl. |n) are the eigenstates and E n the corresponding eigenenergies of the Hamiltonian 
H of the system. The two-point correlation function 



describes correlations between the diagonal matrix elements of a and b in the eigenbasis In Eq. (J2J), (d(E)) AE 

is the mean density of states, given by (d(E)) AE = il/(2irh)f , where / is the number of degrees of freedom of the 
system and ft = J dx5(E — i?(x)) the volume of the constant-energy surface in phase space. The brackets (. . . )ae 
denote a smooth (e.g., Gaussian) energy average over an energy window AE much larger than the mean level spacing 
but classically small, i.e., {d(E)) AE <C AE <C E. The mean density of states determines the Heisenberg time, 

T H = 2nfi(d(E)) AE . (3) 

As shown in Ref. |16| , a second average is required to obtain a self-averaging quantity for the Fourier transform of the 
correlation function R a b(e). We thus introduce the form factor as a time average of the Fourier transform of R a b(e) 
over a time window AT, with AT <C Tjj (for instance, AT = 2%%/AE). Denoting by h(e) the Fourier transform of 
the weight function in the time average, we define 

00 

K ab (T)=(d(E)) AE [ dee-^ T / h h(e)R ab (e). (4) 



This generalized form factor which has been introduced in Ref. ^3] will be the central quantity of this pa- 
per. For definiteness, we consider a uniform average over a time window [T — AT/2,T + AT/2] implying 
h(e) = (eAT/2h)- 1 sm(eAT/2h). 

Setting a = b = 1 in Eqs. J2J) and one recovers the well-known spectral form factor K(T) = Kn(T). The 
correlation function R(e) and its Fourier transform Q have been calculated based on random-matrix assumptions. 
For systems with time-reversal symmetry, the relevant random-matrix ensemble is the Gaussian orthogonal ensemble 
(GOE) and yields the spectral form factor 

K(t) = 2r-rln(l + 2r) , 0<r<l 

= 2t-2t 2 +0(t 3 ) , (5) 
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independent of the dimensionality / of the system. Here, K(t) is expressed in terms of the rescaled time r = T/Th- 
As follows from Snirelman's theorem |17| . the corresponding generalized form factor reads, to leading order in H, 



K ab (T)~a( X )b( X )K(T) (6) 



(see section ITTll below) . Here, the average o(x) of the Weyl symbol a(x) of the quantum observable o is taken with 
respect to the Liouville measure, 

<x) = ^ / dx<5(£ - JT(x)) o(x) , (7) 



see also (|10|) . The phase-space coordinates are denoted by x = (q, p). In many interesting situations a(x) = 6(x) = 
which can always be obtained by shifting a(x) — * a(x) — a(x). This implies, according to Eq. ©, a vanishing K ab (r) 
for fi — ► 0. In this case, our semiclassical methods will enable us to go beyond the result ©. We will show that the 
correction terms to Eq. JHJ are of order 1/Tjj ~ and given by 

oc 

K ab (r) « -±- [2r - 2r 2 + G(r 3 )] | dtCf 6 (t) . (8) 

o 

Here the classical time-correlation function C^ b (t) is defined as 

Ci b (t) = a(x)6«(x t ) with 6 g (x) = &(X) \ b(Tx) , (9) 



where a(x) = 6(x) = is implicitly used, T : (q, p) i— > (q, — p) is the time-reversal map, and x t is the solution of 
the classical equations of motion with initial condition xo = x. The symmetrized form b (x) of 6(x) enters as the 
dynamics is assumed to be invariant with respect to time reversal. 



C. Semiclassical limit 



In the semiclassical approach, the Weyl symbol of the operator a, 



a(q, p) = / dq' > n (q - a q + % ) (10) 
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plays an important role (see, e.g., Ref. [la])- It is a function of the phase-space coordinates x = (q, p) and tends in 
the limit U — > to the corresponding classical observable. In the following we assume that a(x) is a smooth function 
of x. The semiclassical evaluation of Eq. Q for classically chaotic quantum systems yields the generalized Gutzwiller 
trace formula [T^.l20| 



with 



and 



d a (E) = (d a (E)) AE + d° a sc (E), (11) 



n 



(d a (E)) AE = o(x) (12) 




d°a C (E) = ^Rc { > ; Ay exp ( i^L ) }. (13) 

The mean weighted density of states, (d a (E)) * E , depends on the dimensionality / of the system and is determined 
by the average Q of a(x) over the constant-energy surface. The oscillatory contribution, Eq. (|T3|l . is given by a sum 
over classical periodic orbits labeled by 7. The weights A 1 are related to the amplitudes 

w = (T 7 /r 7 )exp(-i7r/x 7 /2) 
V|det(Af 7 ~l)| 
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via the relation A 7 = u; 7 Ay, where T 7 is the period of the orbit 7, r 7 its repetition number, /z 7 its Maslov index, M 7 
its stability matrix, and 

T 7 

A 7 =A(x^,T 7 )ee^- J dto(x?). (15) 

7 

Here, x 7 is the phase-space point on the periodic orbit 7 obtained by solving the classical equations of motion with 
the initial condition Xo , so that x 7 , T = x 7 . The applicability of the semiclassical expression l|13fl to chaotic systems 
with more than two degrees of freedom has been extensively studied in Ref. |2l| . 

Only the oscillating parts of d a (E) and db(E) contribute to the correlation function Substituting Eq. (|T3)l into 
Eq. one obtains 

K ab {T) = J- h£ A^ exp f at ( T - \ {W) 

H \7,7 ^ ^ ^ I AE 

by generalizing the corresponding steps of the semiclassical derivation of the spectral form factor Q . In Eq. I|16|) , the 
delta function with a finite width, 8&t{T'), originates from our choice of time averaging in Eq. It is equal to 
AT -1 if —AT/2 < T' < AT/2 and zero otherwise. The semiclassical formula (|16|l of K a b(T) is the starting point of a 
further semiclassical evaluation. It contains a double sum over terms which strongly fluctuate with energy and poses 
the challenge to approximately compute its energy average. 

An earlier approach to this problem is presented in Ref. [23] where the correlation function (J2J) is considered directly 
instead of the form factor. It was shown that the off-diagonal contributions can be related to the diagonal terms 
yielding the leading order oscillatory term of the corresponding RMT result of R(e) for large e. Furthermore, it was 
pointed out that in the case of a weighted density of states these contributions vanish to leading semiclassical order 
if the microcanonical average of the corresponding observable is zero. In this work we proceed in a different way 
and restrict our considerations to the form factor in the limit of small t = T/Th <C 1, see Subsection IIII Al for our 
corresponding results concerning the weighted density of states. 

In the following section we first discuss the case of the spectral form factor and then generalize our approach to 
include matrix element fluctuations in Section IIIII 



II. SPECTRAL FORM FACTOR FOR /-DIMENSIONAL SYSTEMS 

A. Semiclassical evaluation within the diagonal approximation 

A semiclassical expression for the spectral form factor K(t) = iCn(r) is given by Eq. (jTB)l with Ay — w 7 , By = w 7 
and the rescaled time r = T/Tjj. To leading order in h and r, the double sum over periodic orbits can be evaluated 
by means of the so-called diagonal approximation It is guided by the fact that the contributions from pairs (7, 7) 
of orbits with action differences larger than Ti strongly fluctuate in energy and are therefore suppressed upon energy 
averaging. Hence the main contribution stems from the pairs of orbits with equal actions S 7 — S* 7 . If the system 
has no other than time-reversal symmetry then these pairs are obtained (up to accidental action degeneracies) by 
pairing each orbit 7 with itself or with its time-reversed version 7*. To calculate the corresponding contribution to 
the semiclassical form factor (jltjfl , one has to perform a weighted periodic-orbit average of the type 

<...) p0iT ^i^..> 7 | 2 <5 AT (T-T 7 ). (17) 

7 

This is achieved by means of the following sum rule for periodic orbits in chaotic systems : 

\ T 

d*/(x?)\ ci|dt/(x t )^7M (18) 

po,T 

with T — > 00. On the left-hand side the arbitrary continuous function /(x) is integrated along the periodic orbits 7 
with periods T 7 lying in the interval [T - AT/2, T + AT/2]. The integral on the right-hand side of Eq. JTHJ) is taken 
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FIG. 1: Representation of a periodic orbit 7 (solid line) in phase space with a close encounter together with its time-reverse 
version 7* (dashed line) and its partner orbit 7 P (dotted line). The orbit 7 is characterized by two stretches which are almost 
time reverse of one another. One of these stretches is situated between the two Poincare surface of sections (PSS) perpendicular 
to the orbit shown by the grey squares. The time-reverse map T is the reflection with respect to the plane p — 0. The picture 
should be thought of as a projection of the whole 2/-dimensional phase space on a subspace formed by one momentum and 
two position coordinates. 

along a nonpcriodic "ergodic" trajectory which uniformly an d den sely fills the constant-energy surface. For ergodic 
systems this time average is equal to the phase-space average / (x ) in the large T limit. In the special case /(x) = 1, 
Eq. i|18|) is known as the Hannay-Ozorio de Almeida sum rule . 

For a fixed rescaled time t = T/Tjj, the periods T 7 of the orbits entering Eq. (|16[) diverge as Tjj = 0{ft}~f) for 
% — * 0. This justifies the use of the sum rule l|18|l for the evaluation of the form factor. Since w 1 — iu~< , see Eq. JH}, 
the contribution of the pairs (7, 7) and (7, 7*) to the spectral form factor, neglecting all other orbit pairs, is given 
byi 

K^(t) = 2t((1) po , t ) ae ^2t. (19) 

The factor of two is due to the time-reversal symmetry. It is worthwhile to note that K^(t) agrees with the leading 
term of the RMT result (JSJ for the GOE case; correspondingly the GUE form factor K^ x \t) — r is reproduced for 
systems without time-reversal symmetry. Hence the diagonal approximation explains the universality of the form 
factor for r -C 1 in the semiclassical limit. 



B. Origin of the off— diagonal corrections 

In order to go beyond the diagonal approximation and to explain the agreement of the semiclassical spectral form 
factor with the RMT result (JSJ at higher orders in r, one has to evaluate further terms in the double sum over 
periodic orbits l|lb|) ■ Only pairs of periodic orbits which involve a small action difference of the order of h interfere 
constructively and are not suppressed by the energy average. In a series of papers |a, l3> l3> l2l III] starting from the 
work by Sieber and Richter [g, Q , specific periodic-orbit correlations in systems with time-reversal symmetry have 
been investigated in order to compute the leading off-diagonal corrections to the semiclassical spectral form factor. 
It has been shown that, for hyperbolic dynamics invariant under time reversal, there exists a continuous family of 
pairs (7,7 P ) of periodic orbits with arbitrarily small action differences. These orbit pairs give rise to a contribution 
of K^(t) — — 2t 2 to the spectral form factor. Hence it coincides with the next-to-leading-order term in the RMT 
result . The idea of the approach runs as follows. A periodic orbit 7 is represented by a closed curve in phase 
space. Let us assume that this curve has two stretches which are almost the time-reverse of one another (i.e., they 
are almost identical after applying the time reverse map T : (q, p) 1— > (q, — p) on one of them). In the sequel, we refer 
to such two almost time-reverse stretches of the orbit 7 as a "close encounter" or, more shortly, an "encounter". The 
two pieces of 7 separated by the encounter are called the "left part" and the "right part" . One can associate to 7 a 
partner orbit 7 P by inverting time on, say, the left part, leaving the right part almost unchanged. Hence 7 P follows 
closely the time reverse 7* of 7 on its left part while it follows closely 7 on its right part, as shown in Fig. ^ Such 
a partner orbit 7 P has almost the same action as 7. More precisely, the more symmetric the two orbit stretches of 
7 are with respect to time reversal the closer is 7 P to either 7 or 7* and the smaller the action difference S jP — S T 
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FIG. 2: Configuration space representation of a periodic orbit 7 with a close encounter (solid line) together with its partner 
orbits 7 P (dotted line) for a system with two degrees of freedom. 

Because the periods of the orbits involved in the sum (|16fl are on the scale of the Heisenberg time Th, Eq. ©, one 
expects a large number of encounters on a given orbit 7 in this sum. Thus a large number of partner orbits j p with 
small action differences S 7 v — 5 7 can be associated to any periodic orbit 7. Both, 7 and its associated partner orbit 
7 P share the property to have two almost time-reverse stretches, which are approximately the same for both orbits. 
The partner orbit of j p coincides with the original orbit 7. 

All previous works 0, 0, H, 0, El EI dealing with the contribution of the pairs (7, 7 P ) of partner orbits to the 
semiclassical form factor have been restricted to systems with two degrees of freedom. Then either 7 or j p has one 
(or possibly several) self-intersection(s) in configuration space, which corresponds to the encounter in phase space. 
The right and left parts of the orbit correspond to the two loops formed by this intersection in configuration space, 
see Fig. |2 The right loop is traversed in the same direction while the left loop is traversed with different orientation, 
hence requiring time-reversal symmetry. In order that the two stretches of the orbit near the self-intersection be 
almost symmetric with respect to time reversal, the two corresponding velocities must be almost antiparallel. The 
intersection is then characterized by a small crossing angle s. The orbits 7 and 7 P are distinguished by the fact that 
one has one more self-intersection than the other This is in contrast to the phase-space approach, in which the 
two partner orbits are treated on equal footing. In systems with more than two degrees of freedom, the phase-space 
approach 0, is more appropriate, because for / > 2 the relevant orbits generally do not have self-intersections in 
the /-dimensional configuration space. 

In the following subsections, we will study the spectral form factor of quantum mechanical systems whose classical 
counterparts are Hamiltonian systems with / > 2 degrees of freedom. Furthermore, we consider systems with time- 
reversal symmetry, since only in this case the orbit pairs (7,7 P ) exist. We show that, if the underlying classical 
dynamics is ergodic and hyperbolic, these orbit pairs yield the contribution K^{t) = — 2r 2 to the semiclassical 
spectral form factor, independent of the number of degrees of freedom. Remarkably, the different time scales given by 
the set of Lyapunov exponents {A.;} do not show up in the final result which coincides with the universal second-order 
term of the random-matrix theory prediction JSJ. Our technique strongly relies on the equivalence between the two 
approaches previously developed in Ref. [|| and Refs. [H, |25( to count the number of partner orbits. Therefore we 
present a proof of this equivalence which clarifies the underlying dynamical mechanisms related to the partner orbit 
statistics. Our semiclassical evaluation of the spectral form factor K^ 2 \t) will serve as a basis for the calculation of 
the generalized form factor K^?(r) in Section ITTT1 



Before evaluating the spectral form factor we introduce the notations by very briefly summarizing the necessary 
concepts for dynamical systems [2(|. The classical dynamics of the system is assumed to be ergodic and hyperbolic. It 
maps any phase-space point Xo = x onto the point x t after time t. Hyperbolicity means that all Lyapunov exponents 
are nonzero except the one corresponding to the direction along the flow |26j| . For a given classical trajectory, the 
dynamics in its vicinity can be linearized using the stability matrix M(t, x). The vector Syo(x) = (<5q^, <5pff ) describing 
a small displacement from x perpendicular to the traiectorv|36j within the constant-energy surface is given at a later 
time t by 



This linear approximation is valid as long as ch/t(x) remains sufficiently small. The set of all possible vectors 6yo(x) 
defines a (2/ — 2)-dimensional Poincare surface of section (PSS) at point x perpendicular to the trajectory in phase 



C. Hyperbolic Hamiltonian systems 



Sy t {x) « M(t,x) 5y {x). 



(20) 
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space. The matrix M(i,x) is a linear map from the PSS at x to the PSS at x t . This map is symplectic, i.e., it satisfies 
M T £ M = E, with 



1 

-1 



(21) 



where and 1 refer to the (/ — 1) x (/ — 1) null and identity matrices. Therefore the symplectic product is conserved 
by the dynamics, i.e., <5y t T S 8y[ i=a 5y{f £ <5j/ ' for any two small displacements Syo and Si/q, provided Syt and by[ 
remain sufficiently small. 

The linear stable and unstable directions in the PSS at x are denoted by e/(x) and e"(x). They define vector fields 
which can be found by means of a homological decomposition |2(| of the stability matrix M(t, x). A stretching factor 
Ai(t,x) is associated to each direction. It is defined by 

M(t,x)er< s (x) = A^x)^ eT' s (x 4 ) , (22) 

where the signs + and — correspond to the superscripts u and s, respectively. It is worth noting that Eq. I|22l) is 
not a n eig envalue equation for the matrix _M(t,x), since the vectors e s ' u are evaluated at different positions in phase 
space 37]. In the long-time limit, the stretching factor is related to the Lyapunov exponent A^(x) at point x via the 
relation In | (£, x)| i=s Xi(x)t. On shorter time scales one has to solve the equations of motion 

dAi(£,x) . . , . . . . 

dt =X»(xt)A < (t,x) , (23) 

where Xi( x ) is the local growth rate. In the following, we will assume that the local growth rates are continuously 
varying functions in phase space. In general, Xi( x ) can take negative values in some region of phase space j26j . 
However, by ergodicity, its average Xi( x ) over the constant-energy surface is positive since it is equal to the zth 
positive Lyapunov exponent Aj at almost all points x (i.e., on a set of points x of measure one). The Lyapunov 
exponents of periodic orbits (being of measure zero in phase space) are in general different from the Ai's. They are 
given by 

A?sAiW) = !f*2WM = ^r dtxl(x?) . m 

J 7 J 7 JO 

In hyperbolic systems, the set of vectors {e/(x), e"(x)} spans the whole PSS at x. Hence each displacement vector 
8y(x) can be decomposed into its stable and unstable components, 

/-i 

Sy (x) ee Sy s (x) + Sy u (x) = £ (* (x) e° (x) + u % (x) e V (x)) . (25) 

i=l 

Therefore 5y(x) is determined by the set of stable coordinates {si} and unstable coordinates {ui}. Provided that all 
these coordinates Ui,Si are small enough, the linear approximation l|2UI) can be applied for sufficiently long times t, 
say, up to some time At u 3> A, i _1 (x). By l|22|l and l|25() . the ith unstable component at time t is then equal to its 
value Ui at time t = multiplied by A,(£,x). This leads to an exponential growth of this unstable component as 
\ui\ exp[Aj(x)i] during the time A~ 1 (x) <t< At u . Similar arguments hold for the stable components Sj when going 
backwards in time. This implies an exponential decrease of Si so that the product u, Sj remains constant. For times 
t < At u , it follows from Eq. that 

diii 

-77- = Xt(x t )ttj (26) 

and similarly for Si with \i replaced by — Xi- 

The dynamics uniquely specifies the directions of the e?' u (x.). Due to the symplectic nature of the stability matrix 
M(t, x), they have to fulfill the "orthogonality relations" 

eT(x) T Ee7(x) = e?(xf £e/(x) 

= ^(x) T Se/(x) = (27a) 

for i 7^ j. However, the norms of e^'"(x) can be chosen arbitrarily. In the sequel, we choose these norms in a way that 
their symplectic product gives a classical action S c \ of the system under consideration, e.g., the action corresponding 
to the shortest periodic orbit, 

er(x) T Ee/(x) = S cl %. (27b) 
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encounter region 




FIG. 3: Schematic drawing of the encounter region in phase space. The Poincare surface of section (PSS) is a (2/ — 2)- 
dimensional surface defined at in the 2 /-dimensional phase space by the perpendicular coordinates (5q x , <5p x ). The original 
orbit 7 is represented by the solid line. Also drawn are two segments of the time-reversed orbit Y (dashed lines) which yield 
the two closest intersection points. The corresponding "loop times" t\ are denoted by in and t\2- The displacement vector 
Sy = <5y(x 7 ,ii) points from the original orbit to the intersection points. 



The symmetry of the dynamics with respect to the time-reversal operation T implies M(t, Txt) = TA/(t, x) -1 T, 
where, in the right-hand side, the symbol T refers to the restriction of the time-reversal map to the PSS at x or at 
Txt. It follows that the vectors e^'"(x) can be chosen in such a way that, in addition to Eq. I|27|l. they satisfy 



e7 (Tx) = Tel to , e/ (Tx) = TS? (x) 



(28a) 



and 



Ai(t,Tx t ) = Ai(t,x) 



(28b) 



The relations Q27[) imply that the Jacobian matrix J(x) = d(Sq J -,Sp ± )/d(si,Ui) of the transformation from the 



position/momentum coordinates to the stable/unstable coordinates in the PSS at x is symplectic up to a factor 
i.e., J(x) T S J(x) = — Hence the Jacobian determinant of this transformation equals (— 5 c i)^ _1 . 



D. Encounter region 



As it has been described in Subsection III Bl the pairs of periodic orbits (7,7 P ) which interfere constructively in the 
double sum l|16(l are related to close encounters of 7. Each such encounter involves two orbit stretches of 7 which are 
approximately time reversed with respect to each other. The purpose of this subsection is to give a more precise and 
quantitative definition of the notion of an encounter in the 2/ > 4 dimensional phase space. 

Let us assume that the periodic orbit 7 comes close to its time-reversed version 7* at a point x 7 in phase space so 



that 



Tx 7 



In the following, we choose the time t\ such that Tx 7 +t] lies in the PSS perpendicular to the orbit 



at x/. Thus the small displacement vector between 7 and 7 1 lies in this PSS and can be decomposed in terms of the 
stable and unstable coordinates (see Subsection III C|) . 



fiyfat,t\) 



Tx 7 
1 x t+*i 



/-1 

E 

i=l 



(29) 



( Sl (x 7 ,t 1 )e7(x 7 )+^(x 7 ,t 1 )eTW) 



If one moves from x 7 to x 7 +At along the orbit 7, this displacement vector evolves according to the equations of motion 
and becomes 5y~Atfat > *i) — ^v( x t+At^i ~ 2 At). The displacement 5 y At fat 1 *i) remains small due to the deterministic 
nature of the dynamics if the time At is sufficiently short. In other words, if the two orbits 7 and 7' are close to each 
other at some point in phase space, it takes them a certain finite time until they have significantly deviated from each 
other. 

We define the "encounter region" as the set of all points x 7 +At such that each stable and unstable component of 
the displacement vector 5yAtfat> t\) is smaller than a certain threshold c < 1. The value of c is chosen in such a way 
that SyAtfat ,ti) is given by the linearized equations of motion SijAtfat \h) ~ M(At, x])Syfa] , tj) as long as x 7 +At 
stays within the encounter region, while the linear approximation breaks down outside of it. Therefore c is a purely 
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classical quantity which describes the breakdown of the linear approximation applied to <5yAt(x7, U). As it will turn 
out, the precise value of c is not essential for the calculation of the form factor in the semiclassical limit. This implies 
that a phase-space dependent c(x) does not alter the final result for the form factor. Strictly speaking, c also depends 
on the size of the encounter region, since the corrections to the linear approximation specified above should increase 
with At. However, for smooth dynamics this time dependence turns out to be weak, i.e., logarithmic in \ui\, |sj|, and 
one can show that it does not affect the result for the form factor 

From the definition given above, one concludes that the range of values of At such that xJ +At lies within an encounter 
region is given by — At s < At < At u , where At s , u is defined as follows. Let us denote the ith stable component of 
the time-evolved vector 8yAt(^t,t{) as Si(At;xJ ,t\) and similarly for the unstable components it,. Then At u is such 
that the displacement <5yAt„ is just about to leave the hypercube C — {{si,Ui) : |sj|, < c} meaning that its largest 
unstable component first reaches the value c. A similar definition yields a time At s if going backwards in time, so 
that 

max At s ; x], t\)\} = c and 

max {\ui(At u ;x],ti)\} = c. (30) 
»=i,— »(/-!) 

These implicit equations determine the times At s — At s ({si};xJ) and At u = At u ({ui}; xj ) as functions of the 
components {si,Ui} of the vector Sy(x],t\) defined in and of the point yQ in phase space. The time duration of 
the encounter region, 

t e nc({si,Ui};xJ) = At„({s;}; x^) + At s ({u t };xJ) , (31) 

thus depends on x^ and on all the components {si,Ui}. This time t cnc is clearly invariant within a given encounter 
region. 

The breakdown times At s and At u for the linearization can be estimated in the limit |sj|,|ui| <C 1 by using 
Eq. H30fl and the exponential growth of the unstable and stable components in the forward and backward time 
directions, respectively. With an error much smaller than At u _ s themselves, they diverge like At s w Xj 1 In | 1 1 and 

At u « A^T 1 In Im^T 1 I where j and k are the components for which the maximal values are first reached. For / > 2 
degrees of freedom the presence of the maximum in Eq. (|30|l thus makes the functional dependence of At u _ s and t enc 
on {si,Ui} rather complicated, in contrast to systems with two degrees of freedom. 



E. Partner orbit 



Let us consider an orbit 7 having an encounter at the phase-space location x^ after time t\ as described in the 
previous subsection. For now we assume that the components {sj,itj} of the vector Sy = 6y{j^,ti) are small, i.e., 
\si\, \ v,i\ -C 1. As it will turn out in due course this is the only relevant case for the form factor. We show, by analyzing 
the linearized equations of motion (|20|l around 7 or 7*, that there exists another periodic orbit j p which follows closely 
7 between t and t + 1\ (part 1Z) and follows closely 7* during the rest of the time (part £), i.e., 

7 p J x^ for t < t' < t + 1\ (part K) 

X *' ~ \ Txl t+ti _ t , for t + ti < t' < T 7 + t (part C). ^ 2 > 

Let us denote by Sx-R,i the phase-space displacement between 7 and 7 P at the beginning of part 1Z (time t), see 
Fig-E] To simplify the notations, we do not write explicitly the dependence of the displacement vectors on x^ and t\. 
At the end of part 1Z, i.e., at time t + 1\, the displacement 5xn,i has changed to 5xiz -e . At the beginning and the end 
of part C, the displacement vectors between the time-reversed orbit 7* and the partner orbit j p are denoted by 5x^\i 
and 6x0 e , respectively. Here, C l indicates that one has to invert time on C The vectors 5x are given explicitly by 

Sxn,i = *f ~ , 5x C i,i = x7+ tl - Tx] , (33) 
<5x K , e = xj +ti - x7 +ti and 5x C i >e = xf - Tx] +ti . 

The vectors 8xn,i an d 8x*£i e lie in the PSS defined at xj (see Fig.0J, while 8x-ji,e and 5xc\i are in the PSS at x7 +tl . 
Let R = M(t\,Xt) and L = M{T 1 — t\,xj +t ) be the stability matrices of the parts 1Z and C of 7, respectively. The 
stability matrix of D is given by If = T L~ x T. Since the partner orbit 7 P is assumed to follow closely 7 on part TZ 
and 7* on part C, one can use the linear approximation to evaluate 8x-jz,e and 5xc\ e as functions of Sx-jz^ and Sx^^i, 
respectively, i.e., 



8xn ie = R6xn,i , 8x c , e = L l 5x c ^i- 



(34a) 
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5p x 




FIG. 4: PSS at x t right after part C and before 1Z. The displacement vector Sy points from the orbit 7 to its time-reversed 
version 7*. The deviation of the partner orbit 7 P from the original orbit 7 is described by the vector Sxr^. Shown is only a 
two-dimensional projection of the (2/ — 2)-dimensional PSS. 



These equations determine the two single parts of the partner orbit j p during 1Z and C. In addition, the relations 

&xn,t - Sxc\e = <5y i Sxn, e - 8xa,i = -TSy (34b) 
make sure that the two parts fit together in the encounter region. The set of equations l|34l) can be rewritten to give 
(1 - UR) Sx n ,i = (1 + VT) Sy and (1 - RL l ) 5x C \i = ( T + R) $V • (35) 

Assuming that the determinants of (1 — L l R) and (1 — RL l ) do not vanish, the system of linear equations (|34[l has 
a unique solution. This solution yields the vectors Sx in terms of the displacement Sy. Hence it characterizes the 
geometry of the partner orbit 7 P in terms of deviations from 7 and 7*. 

It is important to note that all points xJ +At within the encounter region lead to the same partner orbit 7 P . This 
means that, when writing Eqs. I|34|) for position xL_ At instead of xjj', the solution is just the vector 8xu,i corresponding 
to xjf shifted along the orbit during time At and similarly for the other vectors Sx in Eq. I|33|l . To see this, let us first 
remark that the time evolution of Sy between t and t + At is determined by the stability matrix M — M(At, xj) via 
SjjAt = MSy. Similar relations hold for the vectors Sx in Eq. I|33|l . The linearization of the equation of motion is, by 
definition, justified within the whole encounter region. The replacement of (x^, ti) by (xJ +At , t\ — 2 At) thus amounts 
to the transformations 

Sy -► MSy 
Sx n .i — > M 5x n .i , 5xn,e 

Sxc.i -> (M t )- 1 5x C ',i , $X£\e 

R -> (M'^RM- 1 , L 
with W = M (At, 

x i+t T -At)- O ne can easily check that the set of equations Q34f) is invariant under these transforma- 
tions. This means that the same partner orbit 7 P is obtained no matter whether x^ or x^ +At was chosen within the 
encounter region. 

Let us first restrict our considerations to the case of long parts 1Z and C. This has to be understood in the sense 
that the linear approximation with respect to the evolution of Sy breaks down at some time between t and t + 1\ and 
similarly, going backward in time, between t and t + t\ — T 7 . This means that t\ > At u and T 7 — t\ > At s . We first 
note that these two conditions actually imply the stronger restriction 

2At u <t l <T 1 - 2At s (37) 

because the displacements Sy at the beginning and the end of parts 7Z and C are related to each other via the time- 
reversal operator T. Formally this can be seen as follows. The displacement Sy^t — fiy(' x -t+At>h ~ 2 At) satisfies 
Syti-At — —TSyAt, as is easily checked with Eq. (|29|) . Let us imagine that At u > t\/2 implying that the linear 
approximation Si/At = MSy is still valid after x^ +At reaches the middle of 1Z. This would imply that |<5j/At| continues 
to increase exponentially with At after time t\/2 until At reaches At u . Such a statement is in contradiction with the 
above-mentioned identity. This shows that At u < t\/2 must hold. A similar argument on part C shows the second 
inequality in Eq. I|37[). 



(M')" 1 Sx^e 

(M>ysx c%e , 

M L M' , 



(36) 
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For a long part 1Z fulfilling l|37(l , the stability matrix R in Eq. Q34a[l is characterized by exponentially large stretching 
factors Aj(x7,ti). Substituting Eq. into Eq. l|33|l and using Eq. J2HJ), one thus finds 

/-i 

5x nti = 5y s = J2 s i ^i( x t) > 

i=l 
/-I 

8xc\i = T6y u = u * e?C rx t') i 

(38) 

<5a? K , e = = - £ *< ^(TxJ) , 

i=l 
/-I 

i=l 

This solution is correct up to first order in the small quantities Si and u;. Terms smaller than and Ui by a factor 
e~ tlX i or e - ^"* 1 - 1 ^ have been also neglected. It means that due to the large lengths of both parts 1Z and C the 
vectors 5x*n,i and T5x^i j describing the partner orbit have to lie very close to the stable and the unstable manifolds at 

x^ , respectively |27j . Furthermore, the points , xj , Tx] +ti , and Tx^_|_ (i form a parallelogram in phase space 0, 0| , 
see Fig. 0] 

It is important to notice that there can be a small set of vectors 8y for which l|37|l does not hold. This is the 
case when either of the parts, say 1Z, is too short jEEJEl- Then the orbit 7 and the time-reversed orbit 7* stay 
close together inside the whole part 1Z so that 1Z is contained within the encounter region. This means that 1Z is 
an almost self-retracing part of trajectory in configuration space. This may happen, for example, in billiards with 
hard walls if one of the reflections is almost perpendicular to the boundary ^JH^]- If there is nopotential or hard 
wall, as in the case of the geodesic flow on a Riemann surface with constant negative curvature (y,0j trajectories 
with almost self-retracing parts cannot exist. If 1Z is contained within the encounter region, the linear approximation 
SyAt = Af(Ai,x7)5y can be applied to 5y — 5y(x],t\) at least up to At = t\. This leads to the additional equation 
RSy — —TSy besides Eq. (|34() . For indeed, following the linearized motion around 7Z, we see that x t and T~x. t +ti are 
interchanged and reverted in time. The solution (|35|l is then Sx-R.i — $y, fi^TZ.e = —T5y and Sx^.i — Sx/h e ~ 0. This 
means that if the time t\ violates the condition l|37|l . the solution of Eq. H34|) does not yield a new partner orbit but 
just the time-reversed orbit 7*. Since the orbit pairs (7,7*) are already accounted for in the diagonal approximation 
(I19H one must only consider intersection points Sy(x];t\) in the PSS which fulfill (|37|l . In other words, the length t\ 
of part 1Z must be large enough so that the linear approximation for 5y\t breaks down for At < t\ and similarly for 
part C. Note that At s and At u are large (of order A J~ 1 In | s ~~ 1 1 and X^ 1 h\\u~^ 1 \, respectively, see Subsection III D|) if 
the components {si, u{\ of 8y are small. 

F. Action difference, orbit weights, and Maslov indices 

The action difference AS* between the orbit 7 and its partner orbit 7 P can be found by expanding the action of 7 P in 
part 1Z in terms of 7 and in part £ in terms of 7*. The derivation of AS* is the same as for systems with f = 2 degrees 
of freedom By using the parallelogram property (|38|l . which is justified since t\ > 2At u and T 7 — t\ > 2At s are 

large, one finds that AS is given in terms of the components {si,Ui} of the displacement 5y(Kj,t\) by 

/-1 /-1 

AS = 5 7 - 5 7 p « 5yT T,5y s = ^ S c \SjUj = y~]Sj. (39) 

i=i i=i 

Thus AS equals the symplectic area of the parallelogram introduced in Subsection III El see Fig. 0] In the last two 
equalities in Eq. H39|) . we have used Eqs. (|27bl I29fl and defined Sj = S c \SjUj. The approximation (|39|l is correct 
up to second order in the small |s<|, -C 1. It is consistent with the concept of the encounter region as it yields 
the same action difference no matter at what position xJ +At within the encounter region it is evaluated. This is 
due to the conservation of the symplectic product under the dynamics. As only small action differences AS* ~ h 
contribute significantly to the semiclassical form factor (|16H , the restriction of the considerations presented above to 
small components |s,|, \ui\ ~ y/Ti/S c \ <C 1 is well justified. 

Besides the two different actions S 1 and S jP entering the semiclassical form factor (|16f) . one must also compare the 
weights w 7 and w 7 p given by Eq. (|14|l . These weights are equal up to small corrections of first order in Ui and s, as 
can be shown in the following way. First of all, for any continuous function /(x) defined in phase space one finds, 
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using Eq. 

| dt'/^V | di'/(x70+ / dt'/(rx2 t+tl _ t ,) (40) 

o t t+ti 

with small corrections of the order of |sj|, \v,i\ ~ yjfi/ S c \. That means that the integral over any function /(x) along 
the partner orbit 7 P is approximately given by integrals along parts of 7 and 7*. The corrections in Eq. H40|) are 
primarily due to the deviations of the partner orbit 7 P from the original orbit 7 or its time-reversed version 7' within 
the encounter region. Obviously, Eq. 140JI yields T 7 ~ T 7 p for /(x) = 1. Similarly we can apply Eq. 140(1 to the 
local growth rates /(x) = Xi( x )j which results into X] ~ X] in view of Eq. (|24() and of the identity x^(Tx) = Xi( x )- 
Hence the Lyapunov exponents of the two partner orbits 7 and 7 P have to be almost equal. Finally, we can also 
identify /(x) with the local change in the winding number of the stable or unstable manifolds which allows for a 
calculation of the Maslov indices [29j. As the winding number of a periodic orbit has to be an integer one finds that, 
for smooth dynamics, the Maslov index of the partner orbit has to be exactly equal to the Maslov index of the original 
orbit ia La, Putting these results together in Eq. I(14jl . one concludes that u> 7 ~ w 7 p. In the spirit 

of a stationary phase approximation we therefore keep only the action difference AS* = S* 7 — S* 7 p in the phase while 
neglecting small differences in the pre-exponential factors in Eq. (|16fl . 



G. Statistics of partner orbits and the spectral form factor 

In the following we show how the orbit pairs (7, 7 P ) specified above determine the next-to-leading-order result for 
the spectral form factor. We assume that the dominant terms beyond the diagonal approximation in Eq. (|16|) are 
due to the systematic action correlations of these orbit pairs. Thus the double sum over periodic orbits can be 
replaced by a single sum over the orbits 7 followed by a sum over all the partner orbits 7 P of 7 while all other terms 
are neglected, i.e., 




where the periodic-orbit average over 7 is given by Eq. (|17fl . All partner orbits 7 P of 7 are characterized by the set of 
action differences {Sj} defined in Eq. i|39|) . Therefore, setting r = T/Th, the sum over the partner orbits in Eq. i|41|) 
can be rewritten as an integral over the Sj's, 

T«-«-(^?u-('g#)) • « 

X -S max (i?) p \ J ~ 1 I 1 AE 

where S max stands for the maximal action difference occurring among the pairs of partner orbits. The density of 
partner orbits 7 P for a given orbit 7 with the set of action differences {Sj} is denoted by d^ _1 A' 7 ({S' :; })/dS'i . . . dSV-i. 
This quantity is the crucial ingredient, and we will show how its periodic-orbit average can be calculated in ergodic 
systems with an arbitrary number of degrees of freedom. In contrast to the case of two-dimensional systems, the 
derivation is significantly more involved because of the higher number of stable and unstable coordinates, Lyapunov 
exponents, and the maximum condition Ij30|l . 

Let us for a moment fix one point x on 7 (to simplify the notation, we omit here the superscript 7 on x and choose 
temporarily the origin of time such that x t = x for t = 0), and consider the PSS V perpendicular to the orbit at x. 
The time-reversed orbit 7* pierces through V many times. Some of these piercings — each of it associated with a 
different time t\ — occur at points Tx tl close to x, see Figs. |3 and 03 Let p 7 ({si, Ui}; x) d?~ 1 s d^~ 1 u denote the number 
of such intersection points, with stable and unstable components of Sy = Tx t] — x lying in intervals (sj,Sj + dsj) 
and (ui,Ui + dui), respectively. We exclude from p 7 all points Tx tl violating the condition 2At u < t\ < T 1 — 2At s 
since they either do not exist at all or do not give rise to a distinct partner orbit. Thus we have the density of valid 
intersection points, 

T 7 -2At s . 1 

p 7 ({si,Ui};x) = / dtiS ([x- Tx tl ]||) JJ S ([x - Tx tl ] u ,i - Ui)S([x- Tx tl ] Sji - Si) . (43) 
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FIG. 5: Sketch of the PSS as it is shifted along the periodic orbit 7 (solid line). Three pieces of the time-reversed orbit 7* 
are represented by dashed lines. If the PSS moves with the flow in phase space from to xj 2 all intersection points of the 
PSS with 7 1 change their positions according to the linearized equations of motion 11201 . Note that not only 7* but also 7 itself 
could come close to at a later time. However, we have not include this in the sketch above. 

Here, the first delta function ensures that Tx tl lies in V, i.e., that the coordinate [x — Txtjii of x — Tx t] in the 
direction parallel to the orbit vanishes. The lower indices s and u indicate the stable and unstable components of the 
vectors inside the square brackets. 

In order to determine how many partner orbits 7 P of a given fixed orbit 7 exist with a given set of action differences 
{Si}, one has to count the number of corresponding encounter regions of 7. As explained in the beginning of this 
section each of these encounter regions can be associated to a displacement vector Sy or to its corresponding time- 
evolved <5j/At, —At s < At < At u . Therefore we consider the dynamics within the PSS Vt at x^. It can be parametrized 
by means of the stable and unstable coordinates of the different vectors 5y(x] ,t\) associated to different times t\. As 
the PSS is shifted following the phase-space flow along the orbit 7, the stable and unstable coordinates of each such 
vector change leaving only the products Sj = S c \ SjUj invariant. The vector SyAt = ^2/At(xt, t\) S Vt+At corresponding 
to a fixed t\ thus moves, as Vt+At is shifted by increasing At, on a hyperbola as long as x^ +At remains within the 
encounter region, i.e., for —At s < At < At u , see Fig. EI Since the number of partner orbits is equal to the number of 
encounter regions one has now to count each encounter region exactly once. This can be achieved in two alternative 
ways: 

1. One can measure the flux of vectors Sy through the hypersurface defining the end of the encounter region (see 
Fig. [SJ). According to the definition of the encounter region given in Subsection III Dl this hypersurface consists 
of the faces dCf = {(si, Ui) : \si\ < c, \ui\ < c, Uj = ±c}, j = 1, ...,/— 1, of the hypercube C = {(si, Ui) : \si\ < 
c, \ui\ < c}. The union of all these faces defines a (2/ — 3)-dimensional closed hypersurface dC contained in the 
(2/ — 2)-dimensional PSS. The corresponding flux is obtained by multiplying the density p 1 with the component 
iij of the velocity of the vector Sy in the direction normal to dC^ 1 . This velocity is given by iij = Xj(^-t) c ; see 
Eq. l|2rj|l . Integrating along the orbit we obtain 



The last product of delta functions restricts the action differences to the values {Si}. It should be noted that, 
since the local growth rate Xji^t) can take negative values for some times t, the vector Sy may also re-enter into 
the hypercube C through some face dCf with a negative normal velocity iik — Xfc( x 7) c (with possibly fc 7^ j). 
However, since Sy increases exponentially with time at large times, there is one more passing of Sy through dC 
in the outwards direction (ii > 0) than in the inwards direction (it < 0). The contributions of all subsequent 
passings then mutually cancel each other in Eq. (|44|l . Hence, for each encounter region, only the first crossing 
of dC at time At = At u is accounted for, as required. Let us also mention that if we had taken any other closed 
hypersurface contained in the hypercube C instead of dC the same result would have been obtained. This is 
because the dynamics conserves the number of points in phase space and thus the number of vectors Sy. 

2. An alternative version of Eq. I|44|) . treating all points within the encounter region on equal footing, can be found 
as follows. Every vector Sy is counted as long as it remains within the hypercube C. Therefore one has to include 




T- 




(44) 
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FIG. 6: Schematic drawing of a projection of the Poincare surface of section (PSS). The flow of intersection points (black filled 
circles) is represented by the thin arrows. For long enough time the unstable component Uj grows while the stable component 
Sj shrinks with Sj = S c \SjUj being constant. There are two ways to count the intersection points. Either the flux through the 
Uj — c surface (dotted line) is considered, as in Eq. 11441 . or one counts the number of points in the volume of the encounter 
region (dashed area) and normalizes that by the time each point spends in there, as in Eq. 1451 . 



the additional factor of l/t cnc , since per definition l|31|) . i cnc is approximately the time each vector 6y spends 
within that hypercube. The density of partner orbits 144|) can thus be rewritten as 



dSi...dS/_i J J tcnc({Si,Ui};X t ') V AJ - 

-c ' 

More precisely, this expression can be derived as follows (also see the Appendix J . Consider first only the 
contribution of the encounters at after time t\ — 2t, for a fixed t\ and an arbitrary time t. The time duration 
t onc of the encounter and the product itjSj of the stable and unstable components of the vector 5y(xJ,t\ — 2t) are 
independent oft as long as xj stays within an encounter region, i.e., while Si and Ui vary within (— c, c). The 
time spent by the point characterized by 6y(~x.] , t\ — 2t) within the hypercube C is approximately equal to t cnc 
in the limit |it,|, \si\ ~ yJh/S c \ -C 1, where t cnc is large (of order ln?i). Indeed, although possible re-entrances of 
Sy into C (see above) may increase the total time spent by Sy inside C to a value greater than t cnc , the relative 
error made by approximating it by t cnc is small. Using also the fact that the above-specified encounter regions 
are disjoint (if they were overlapping they would define one bigger region), it follows that the r.h.s. of Eq. I|45|) 
gives approximately the density of those encounter regions with respect to the action differences {Si}. The 
result l|4*5|) then follows by integrating t\ over all possible values. 

Actually, expressions (|44ll and H45[) are equivalent. Using the fact that the number of vectors Sy is conserved by 
the dynamics, one can transform the integrals over the hypersurface dC into an integral over the entire volume of the 
hypercube. More details on this proof of the equality between Eq. (|44|) and Eq. (|45|l are given in the Appendix B It is 
important for what follows to note that this equality still holds if the range of integration of t\ in Eq. (|43|) is replaced 
by the larger interval (0,T 7 ) which corresponds to the additional inclusion of intersection points Sy that cannot be 
associated with a partner orbit. The reason for the introduction of the two different expressions for the number of 
partner orbits is because of its crucial importance from a technical point of view. We will apply either Eq. 144(1 or 
Eq. i|45|) depending on which one can be calculated easier. It will turn out that this allows for a major simplification 
of the derivations to follow. In particular, the complicated analytic structure of t enc ({si, Ui}; x), see Eqs. (|3U|I and 
(loT^ . docs not directly enter the calculations. 

The periodic-orbit average of the density of partner orbits in Eq. (|44|l or Eq. I|45|) can be transformed into an 
average over the constant-energy surface by means of the sum rule l|18|) . After this step has been performed, the 
density p 7 has to be evaluated at arbitrary points x on a set of measure one inside the constant-energy surface, instead 
of taking the points x^ belonging to periodic orbits as the arguments. For such points x one can neglect classical 
correlations between x and Tx tl for 2At u < t\ < T — 2At s . This is because A" 1 At u ^ s <C T in the relevant limit 
\si\, \ui\ ~ y/h/S c \ < 1, T~ Tjj ~ Ti 1 ~f . More precisely, ergodicity allows one to approximate the time integral in 
Eq. il4*5|) by a phase-space average, 
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IS \ \ ~ r-2At fl ({ai};x)-2At M ({Mi};x) f , , „, , m/[ , 
pds^-ujjx) w / dx d(E - H(x))5 {[x-Tk\\\) 

/-i 

J] 5 ([x - Tx'] u>i - ui) <S ([x - Tx'] Sii - Si ) . (46) 
»=1 

Since the Jacobian of the transformation (Sq^-fSp- 1 ) — > (sj,Uj) gives a factor S^ -1 this yields 

p({ Sl , Ul };x) « p lcad + p COT1 '({ Sl , Ul };x) (47) 
r-^2i cnc ({ S „ Ml };x) . 





Therefore p is given by a leading contribution plus a small correction term due to the exclusion of short times U violating 
condition 1)3 7|1 . The corrections to the ergodic approximation are not written in Eqs. (|46J) and 1)47(1. Although they 
may be bigger than p corr , one expects them to be strongly reduced after averaging x over the constant-energy surface, 
as required by the sum rule. This is not the case for the correction term p corr , which, as we shall see now, determines 
the form factor. 

Indeed, if only the leading term p lead in the density (|47l) is considered, one finds that the form factor l|42[) vanishes 
in the semiclassical limit for the following reason. As p lead does not depend on x and {si, Ui}, its contribution to the 
density of partner orbits can be most easily calculated by means of Eq. . It yields 



Here we have used the identity Xj( x ) = ^j- If this result (|48(l is inserted into the expression for the form factor l)42)l . 
one obtains K ^ (r) = due to the energy average. 

Therefore the small correction term p corr given in Eq. (|47|) is of crucial importance. To determine its contribution 
to the form factor, it turns out to be technically favorable to use expression l|45[) instead of Eq. 1441) for the density of 
partner orbits. The reason is that the two appearances in Eqs. (|45() and l|47|) of t cnc mutually cancel. Inserting p corr 
from Eq. ljT7|) into Eq. (gSJl, one finds 

/«)\ H - _ 2 / 1 rf m (^) (49) 
\d5 1 ...d5 / _ 1 / poT - n l\ ln { |4| ) ■ (49) 

The result © together with Eq. |@SJ| gives the correct asymptotic form of the averaged density of partners in the 
limit h — > 0, r = T/Th fixed. Since the leading term l)48|l gives a vanishing contribution to the form factor (|42|l . only 
the correction H49|) determines the final result, 



K^r) c -2^ TJ J- S / d* e„ P (if) 1„ (0) j c . (50, 

This result is universal, i.e., it does not contain any information about the set of Lyapunov exponents {Xi} or 
the constant c defining the encounter region. Thus as our first major result we find that the next-to-leading order 
correction beyond the diagonal approximation agrees with the BGS conjecture independently of the number of degrees 
of freedom the system possesses. 



III. MATRIX ELEMENTS FLUCTUATIONS 



The aim of this section is to evaluate the generalized form factor K a b(T) defined in Eq. J3J based on the method 
developed in the previous section. This form factor describes the correlations of the diagonal matrix elements (n|a|n) 
and (m\b\m), corresponding to distinct energies E n and E m , of two given quantum observables a and b. We assume 
that a and b have well-behaved classical limits given by smooth Weyl symbols a(x) and 6(x). 
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A. Leading term 

To zeroth order in Ti, the form factor is given by 



K a6 (T)«a(x) b(x)K(X) . (51) 

Actually, Snirelman's theorem for classical ergodic flows implies that the Wigner functions of almost all eigenstates 
\n) with energies E n converging to E are uniformly distributed over the constant-energy surface -ff (x) = E in the 
semiclassical limit. Equivalently, this means that the matrix elements (n|a|n) converge to the average |J7||. 



(rt|a|rt) — > a(x) 

as % — ► , n — ► oo such that E n K, E . (52) 

It is worthwhile to mention that this is only true for eigenstates |n) of the quantum Hamiltonian pertaining to a 
"set of density one" J^]. Heller's scars [^3 are prominent examples of "exceptional" eigenstates violating Eq. 152JI . 
Choosing, e.g., a Gaussian weight of width AE in the energy average in Eq. @, one can express K a i,(T) for T > AT/2 
as 

K °»w = jink — ^^fE^i^^^i^)^^^" 1 ^"^^^™-^) 

/ (E n +E m -2E) 2 \ 
X6X H 8AE 2 J' (53) 

Since all functions of E n and E m vary noticeably on the scale AE (<i(i?))^£, the eigenstates not belonging to 
the "set of density one", such as scars, have a negligible contribution to the sum in Eq. I|53|l . One can then replace 
(n\a\n) (m\b\m) by the product a(x) 6(x) and move this factor out of the sum. This yields Eq. I|51|l . which is therefore 
a direct consequence of Snirelman's theorem. 

To obtain information on matrix element fluctuations, one thus needs to study the semiclassical corrections (next 
term in power of h) to the leading behavior (|51|l of the form factor. Let us define 



a=a-a(x)l , b' = b-6(x)l, (54) 

so that the associated Weyl symbols a'(x) and 6'(x) average to zero. Then the form factor K a b(T) is related to 
K a 'b' (T) by the formula 

K ab (T) - a (x) 6(x) K(T) = (55) 
K a , v {T) +ti^)K w (T) +W)K a >i{T) • 

Comparing with Eq. (|51|l . one sees that the r.h.s. of Eq. (|55fl vanishes as Ti — > 0. The purpose of the two next 
subsections is to estimate the first term of the r.h.s., which turns out to be proportional to 1/Th = 0{fi^ 1 ). We 
start with the diagonal contribution of pairs of identical orbits (modulo time-reversal) to K a n,'{T) and then include 
the pairs of correlated orbits (7,7 P ) studied in Section [H] We restrict our derivation to the case of observables with 
vanishing mean, i.e., a(x) = 6(x) = so that a = a' and b = b' . Therefore we shall not be concerned further in this 
paper with the second and third terms in Eq. (|55J) . 



B. Correction term within the diagonal approximation 

Let us first consider the semiclassical correction to the leading term (|B 1|) within the diagonal approximation. This 
correction has been already studied in Refs. 0, 0, Q, 0|. However, we will argue below that the results of 
Refs. 0,0,0,0 can only be applied to observables a(x) or 6(x) independent of the momentum p. We treat here 
the more general case of smooth observables a and b depending on both the position q and the momentum p, by 
following the lines of Subsection IIC of Ref. [13] . 

Retaining only the contribution of those pairs obtained by pairing each orbit with itself or with its time-reversed 
version in the double sum l|16(l . the semiclassical form factor can be written as 
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K a^T) = ±- ( f'" * a(x?) ( P dt" 6(x7„) + d*" b(Tx]„)) ) . (56) 




po,T 




Substituting t' = t" — t and using the periodicity of 7 yields 

r ^ dt' \ 

^a(x7)6 s (x7 +t ,)) (57) 

7 ' po,T 

with 6 s (x) as given in Eq. 0. 

We now assume that the classical dynamics is sufficiently chaotic so that the time-correlation function of the 
classical observables a(x) and b s (x) decays faster than 1/t to zero. In strongly chaotic systems all classical correlation 
functions of smooth observables decay exponentially, as a result of a gap in the spectrum of the resonances of the 
Frobenius-Perron operator (all resonances but the one corresponding to the Liouville measure are contained inside a 
circle of radius strictly smaller than unity) [2(j. The mixing property makes sure that the time-correlation function 
C^ b (t) tends to zero in the large-i limit, but is still not strong enough for our purpose: it does not imply that C^ b (t) 
can be integrated from to 00. 

Applying the sum rule (fT%f to Eq. ifBTjl gives 



K™(T)~— dt'C s ab (t<) (58) 



1 h J a 

in the limit h — > with r = T/Th fixed. If a(x) or 6(x) is a function of the position q only, then a s (x) = a(x) or 
& s (x) = 6(x), respectively. As a result, C^ b (t) = C a b(t). In such a case, Eq. I|58|l coincides with the result of Refs. |l^ 
and 13]. In the opposite case, the mean C^ b (t) of the correlation function of a(x) and 6(x) and the correlation 
function of a(x) and b(Tx) must be considered. 

As noted in Ref. |l3j . some chaotic systems which fulfill the mixing property, such as the symmetric stadium billiard, 
exhibit algebraic decays of correlations in 1/t. Then the integral in Eq. Ij58(l diverges and the form factor K^ b is of 
order T^ 1 hxT H instead of T^ 1 . 

Let us also mention that, for chaotic dynamics, the integrals J Q T dta(x t ) and J Q T dtb s (x t ) of the observables a(x) 
and b s (x) along pieces of (nonperiodic) trajectories of time T, thought as function of the initial point x, may be 
often considered as Gaussian random variables with respect to the Liouville measure for large T's [26|. These random 
variables have a system-specific covariance 2T J* °° dtC^ b (t), which is thus also related to the fluctuations of the 
diagonal matrix elements of a and b as given by K ab (T). 



C. Contribution of the partner orbits 

The contribution of the partner orbits to the semiclassical form factor l|16|) is obtained by inserting the product 
A 7 B-yP of the integrals of the classical observable a(x) along 7 and of the observable b(x) along 7 P in front of the 
exponential in Eq. I|41|) . The forthcoming calculation is simplified by noting that if j p is a partner orbit of 7, then its 
time-reversed version 7 PjI is also a partner orbit of 7, with the same action. This is because if one exchanges the role 
of parts 1Z and C in the definition l|32|) . the corresponding partner orbit is just 7 P ' 4 . As a result, one may equivalently 
insert A 1 {B 1 p + B^ P .i)/2 in front of the exponential in Eq. I|41|). instead of A 1 B lV . The mean (B^ P + B jP ,i)/2 is the 
integral of the symmetrized observable b s (x) given by Eq. © along 7 P . It can be estimated by applying Eq. I|40() to 
b s (x) and using b s (x) = b s (Tx) together with the periodicity of 7. This yields 




and reflects the fact that the two orbits 7 and explore almost the same phase-space regions as the two partner 
orbits 7 P and 7 P, \ Hence the generalization of Eq. (|42|l reads 
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By using Eq. I|44|l and substituting t' — > t'" = t' — t before applying the sum rule JTSJl, one finds that the leading 
contribution p lcad to the density in Eq. 1(47(1 yields 



' d'-^my) ^ dt> f-< df a(x7)5 5 fx7 A (load) 
ds 1 ...ds f .J Q tJ T , a ^ b (x ^"7 

\ I po,T 

f-l /"I/ 



E II ln (^f ) [ j\t" X] {K)a(* t ,,,)bS{* t ,,, +t , l ) . (61) 



ft 

3=1 i^j 



Employing ergodicity, the integral over t'" can be approximated by a phase-space average and yields 
Txj(x) / °°di"C*f fc (i"). Thus inserting Eq. (J^TJ into Eq. ((HDJ gives ^(r) = 0. As in Section EH the contribu- 
tion to the form factor of ( o load thus vanishes. By using Eq. ((45(1 . we obtain the contribution of the small correction 
term p( corr ) in Eq. gJJ, 



/ ^-^({gj}) / T " ^ f T ^ a(x , )W A corr 

\ / po,J 

»^'^II'»(^f) [ <*"[ d t »„(*,»»*(x,. + ,) (62) 

since the dependence on i enc ({si, i^}; x) in Eq. H45|) and Eq. (|47|l mutually cancels. The average a(x t /»)& s (x t »/ +t ») 
equals the correlation function C^ b (t"). It follows from Eq. 150|) that, as h — > 0, 

K^(t)^-2t — J o dtC s ab {t). (63) 

Remarkably, one obtains for the leading off-diagonal contribution the same result as for the diagonal approximation, 
with 2 replaced by —It as in the spectral form factor. In particular this means that the classical correlations enter in 
exactly the same way via the correlation function C^ b (t). Assuming that only identical orbits modulo time-reversal 
symmetry and pairs of partner orbits (7,7 P ) contribute to the semiclassical form factor 1(16(1 up to order t 2 included, 
this yields 

K ab (T) « ^r(K(T) + G(r 3 )) J~ diC&(t) (64) 

as announced in the introduction. This result holds if the correlation function C^ b {t) decays faster than 1/t as t — > oo, 
in order that the upper integration limit T = tTh may be replaced by oo. It is valid for observables a and b such 
that a(x) = 6(x) = only. 

If a(x) ^ or b(x) ^ 0, Eq. 1(55(1 must be used and the second and third terms in the r.h.s. of this equation have 
to be estimated. Repeating the above calculation for these terms, one finds that both vanish in zeroth order in h, 
thus being in accordance with Snirelman's theorem. For instance, within the diagonal approximation, Eq. I|56|) gives 

Kf£\(T) — 2T II 1 (J 7 dt a'(x7)) po ,T ~ 2r a'(x), which is zero since a'(x) = 0. The leading contribution in h is thus 
governed by the finite-time corrections to the sum rule ((18(1 . Similarly, replacing a by a' and b by 1 in Eq. 1(62(1 . the 

(2) 

second integral in the second member becomes Ta'(x), which means that KX{(T) fa up to higher-order corrections 
in the sum rule. One concludes that our method does not allow us to estimate K a >i as % — > beyond the leading 
order in ti. For systems with exponential decay of classical correlation functions, it is not irreasonable to expect that 
the finite-time corrections to the sum rule (|18f) are exponentially small in T. In such a case K a i\ and Kw would be 
negligible with respect to K a n,>, which is of order h^ 1 by Eq. 1)64(1. 
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IV. SUMMARY AND OUTLOOK 



In this work we presented a semiclassical evaluation of the generalized form factor K a b(r) going beyond the diagonal 
approximation. We first considered the spectral form factor K(t) — Kh(t) for systems with more than two degrees 
of freedom, i.e., / > 2. We proved that the leading contribution due to pairs of periodic orbits with correlated actions 
is independent of / in agreement with the RMT prediction. An important step in our calculation was to show the 
equivalence between the two different approaches for counting partner orbits which were independently developed 
in Ref. Q and Refs. 0, for two-dimensional systems. Based on these results for the spectral form factor we 
then investigated the generalized form factor K a b{r). In this case we were able to show a universal dependence of 
K a b{r) on the rescaled time r. Furthermore, we found that the contribution of the partner orbits depends on the 
classical time-correlation function C^ b (t) in exactly the same way as in the diagonal approximation, see Eq. Ij64(l . An 
interesting open question is to prove (or disprove) that this is still the case at higher orders in t or even for arbitrary 
large r > 1. In such a case one could get rid of the error term C(r 3 ) added to K(t) in Eq. (|64fl . 

Our semiclassical treatment of the generalized form factor beyond the diagonal approximation can in principle 
be extended to other physical observables containing matrix elements in chaotic systems. This includes expressions 
where transition matrix elements play a role [3lll33| (e.g., dipole excitations in quantum dots and linear response 
functions for mesoscopic systems |34j with applications to transport, magnetism, or optical response. So far, nearly 
all semiclassical approaches to such quantities have been relying on the diagonal approximation, as long as additional 
averages are involved. A notable exception is the calculation of the weak localization correction to the conductance 
in Ref. showing an important contribution of the partner orbits. It would thus be of great interest to study the 
corrections to the diagonal approximation in the various response functions appearing in mesoscopic physics. 
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APPENDIX: TRANSFORMATION OF THE SURFACE INTEGRAL INTO A VOLUME INTEGRAL 

In this appendix we prove the equality of the two different approaches for counting the partner orbits based on 
Eq. 144|) and Eq. I|45|) . respectively. To this end we show that an equality of the general structure 

T T 

Jto Jwn$rA = J dt J dAsy g{5y, t) 8y(5y, t) (A.l) 

V ' dV out 

holds under the conditions which are relevant for the statistics of the number of partners. Here, Sy is a vector 
in a multidimensional space, e.g., the (2/ — 2)-dimensional PSS. The volume element in this space is given by 
dV^ = d^ _1 ud^ _1 s while dAgy- characterizes the surface element. The left hand side of Eq. (|A.1|I thus contains 
an integral over any (2/ — 2)-dimensional volume V in the PSS. Inside V we follow the time evolution of a density 
field g(Sy,t); the corresponding velocity field is denoted by Sy(Sy,t). As Eq. (|A.1(1 is applied to the PSS following a 
periodic orbit of length T, we can assume periodicity such that g(Sy,t) — g(5y,t + T) and Sy(5y,t) — Sy(Sy,t + T). 
Due to current conservation the density is constant along the flow, i.e., g(Sy, 0) = g{Sy t , t) or g(5y t ,t) — 0. The time 
ty(Sy,t) in Eq. (jA.lfl is defined as the total time a point spends in the volume V if it starts at time t at position 8y 
and moves until time t + T. If the volume V is chosen to coincide with the hypercube C defining the encounter region, 
see Subsection III Dl then t\> is approximately equal to the time t enc , Eq. (|31fl . The surface of V is decomposed as 
dV = <9Vin + <9V out . Here, <9V in / out stands for that part of the total surface through which the flux defined by g and 

8y enters or leaves V in the long-time limit, respectively. More precisely speaking, the total flux between time and 
T through any piece of <9V ou t must be positive. 

For the proof of relation (|A.1|I let us first consider the case where the total density g(Sy, t) is given by a single point 
starting at Syo, i.e., gi(Sy,t) — 6(5y — Syt). Then the time iy is given as 



t+T T 

" Sy ' 



tv{5yut) = J dt' Q v (Sy t ,) -Jdt' J dV sp gi(Sy',t') 
t o v 

= t v (6y ,0) , (A.2) 
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where Q\>(5y) equals 1 if Sy € V and zero otherwise. In Eq. (|A.2|I we made use of the periodicity of the motion. We 
then obtain for the left hand side of Eq. (|A.1() 



T 



/ 



* /d%f#$= (A.3) 



o v 



tv(8y,t) 

T 

^-^ Jdt J d% Qi(5y,t) = 1 . 



t v (6y 

In close analogy we thus find that if the single point density is replaced by g{Sy, t) = J2i Qi(Sy, t) which represents an 
arbitrary number n of points given by their initial conditions then the left hand side of Eq. ijA.lfl just gives the total 
number of particles n that pass V during one period. But this is exactly what the right hand side of Eq. IjA.ll) gives. 
It just measures the outgoing flux through the surface of V between time and T which also yields the total number 
of particles n because the particle number is conserved. 

Finally we also note that the density g(5y, t) is not restricted to a sum of 6 functions. Each of these 5 functions can 
also be multiplied with any function g(Sy, t) that is constant when following the flow within V, i.e., g(8y*o, 0) = g(Syt, i). 
In the context of Subsection lll (Jl g could, for example, be any function of the action difference as in Eqs. (|44|l and 145|l . 
In this case the density g entering Eq. I|A.1|) can be considered as a weighted density g — pg. 

If all local unstable growth rates Xfc( x ) are non- negative one can directly identify iy = t onc and thus the equality 
(|A.ip means that Eq. 14411 exactly equals Eq. H45|) . On the other hand, if these local unstable growth rates assume 
negative values in certain areas of the phase space then this implies that the unstable components of a displacement 
vector can also decrease on short time scales. This would lead to a multiple entry of the same point into the 'encounter 
region' characterized by V — C. In this case the relation (|A. 1|) means that Eq. 145|) is asymptotically equal to Eq. 144|) 
as the length t cnc becomes large so that |ty — t onc | <C t cnc or similarly £y ~ t cnc . 
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